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The theory for the diffraction of spherical electromagnetic waves by a finitely conducting 
spherical earth was developed from Maxwell's equations by Watson [1918] and the intricate 
computation details were later worked out by van der Pol and Bremmer [1936] as the now 
classical series of residues. Two aspects of this computation present considerable difficulty, 
especially at low frequencies : 

(1) The calculation of the height-gain factor which takes account of an elevated trans- 
mitter and/or receiver. 

(2) The evaluation of the special roots, t=Ts, of Riccati's differential equation, 

g-2.^. + l=0, 



near the circle of convergence, l52r| = J^. 

These analytic difficulties are avoided with the aid of modern analysis techniques 
applied to a large scale electronic computer. Hankel functions of the first and second kind 
of order one-third and two-thirds are calculated by numerical integral methods and then 
used with iteration to solve Riccati's differential equation. The amplitude and phase of 
the spherical radio wave diffracted in the vicinity of the earth with various altitudes above 
the surface of the earth, of both the transmitter and the receiver, are then calculated by a 
summation of the series of residues. 



1. Introduction 

The field of the diflfracted spherical wave in the 
vicinity of a finitely conducting spherical earth has 
been discussed in detail as the ground wave by 
Watson [1918], van der Pol and Bremmer [1937], 
and the particular case of a vertically polarized 
Hertzian dipole source, Er, volts/meter; i.e., the 
electric field directed radially with respect to the 
center of the earth, can be written [Johler, Kellar, 
and Walters, 1956], 



Er = 2E,rFr, 



(1) 



where, 



E,r-^-^Iolexp{i[k,d-o^t]}, (2) 



where, d is the distance, meters, /=co/27r the fre- 
quency, cycles/second, IqI is the dipole current 
moment, ampere-meters, exp {—ioit) implies a wave 
varying harmonically in time, ^, in a medium with 

wave number, ki=- r/i, where c is the speed of hght, 

c^S (10^) m/sec, 771 is the index of refraction of air, 
Tyi'^l. The attenuation function or secondary 
factor, Frj which takes account of the earth's 
curvature, finite conductivity, and dielectric con- 



1 Sponsored by Air Force Cambridge Research Laboratories, CSO&A 58-40 
(3/31/58). 

2 An earlier version of this paper was presented at the 1961 Spring meeting of 
the International Scientific Radio Union (URSI), Washmgton, D.C. (May 1961) 



stant, can be evaluated as a summation of the 
classical series of residues, 



L ^J S = o 1 I 



X 



»p{'[»'«''-'5+i;+i]}' 0) 

where a is the radius of the earth and a is a factor 
which takes account of the vertical lapse of the 
permittivity of the earth's atmosphere. 

The parameter, b=be for vertical polarization is 
defined, 



5.=- 



'K2 1 



(kio) 



•[i-'J 



(4) 



€' 



wit\ik2=-\ e2+i — ~ r radians/meter, €2 the dielec- 



tric constant relative to a vacuum, a the conductivity 
of the medium, mhos/meter, /Xo the permeability of 
free space, henrys/meter, and t=Ts describes the 
special roots of Riccati's differential equation, 



d8 



-25V + 1=0, 



(5) 
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which define the positions of the poles at which the 
residues for the series 8=0, 1,2, 3 . . . are calculated. 
In the case of horizontal polarization, E^ may be 
replaced by Hr (vertical magnetic field) providing 
6e is replaced by bm where, 



C' 1 ? 

/t2 



(6) 



The isolated factors, JsiK), .fs{h2), which are in 
each term of the series describe the effect of elevating 
either the transmitter, hi, or the receiver, ho, or 
both, and can be evaluated without further ap- 
proximations [Johler et al., 1956] as follows: 



m)= 



iha) 



2/3 



2ha"^ 



-2t, 



1/2 






(7) 



It is often more convenient to use some sort of 
modified Hankel function, hi(z) [Furry, 1945], for 
which values of complex argument (z=x-[-iy) have 
been tabulated. 



h 



fsihy- 



■{(if [(..„=-. g^-.,.]} 



(8) 



hi{-(2y^'r,} 



The roots, r^, are found [Johler, Walters, Lilley, 
1959] by expanding the Riccati differential equation 
in a power series, 






\S'rKh 



-'^s,0 



^ 1 r d^rs 



5V|>|, 



87", 



(9) 



(10) 



where the hmiting roots, Ts,o and Ts,a, for 8e=0 and 
5g=z 00 have been tabulated [Johler, Walters, Lilley, 
1959]. 

The power series in 8"^ and 6~^ both seem to lose 
precision near the circle of convergence, \d^T\=^, for 
a finite number of terms. Although a considerable 
number of terms have been determined [Johler, 
KeUar, Walters, 1956] the calculation of more terms 
becomes very laborious. Furthermore, the tabula- 
tion of the Hankel functions [Furry, 1945] does not 
seem to provide an adequate range of values, es- 
pecially at low frequencies and indeed is not especially 
suitable for application to large scale computers. 

Since the Riccati differential equation can be 
expressed as a ratio of Hankel functions, 



m]i\ 



4(-2r,) 



3/21 



■^^=27.m%m-2ur 



exp 



[^f]+^^=0' 



(11) 



a basic technique to evaluate Hankel functions 
^j/3?2/3 (^) would permit an iterative solution of 
Riccati's differential equation and at the same time 
could be employed to evaluate the Hankel functions 
H[]l (z) in the height gain factor, J,(h), The 
development of such a technique is the prime task 
of this paper. 

2. Evaluation of Hankel Functions by 
Gaussian Quadrature 

Bessel functions of order v are solutions of the 
differential equation. 



^^g+^t+(^-^^)^=o 



dz 



dz 



(12) 



where v and z are real or complex. The Bessel func- 
tions of the third kind, known as Hankel functions, 
are linearly independent solutions of eq (12) and are 
denoted by H^^^ {z) and Hl^'^ {z) where v=\ or ^ 
for the purpose of this discussion. Hankel functions 
have complex values for a real argument, but are 
real for r+^ H,^'' (iy) and i-^^+^> m'\ {-iy) when 
y is positive. These functions are important in 
applications because they are the only Bessel func- 
tions that vanish for an infinite complex argument, 
H\}^ {z) if the imaginary part is positive and H^^^ 
{z) if it is negative [Jahnke and Emde, 1945]. 

Integral representations for Hl^^ {z) and Hl^^ [z) 
are as follows [Jeffreys and Jeffreys^ 1956]:^ 

where the expression of the limits in eq (13) means 
that the path is taken from to — od by way of i. 

Integrals with real limits are derived by selecting 
a path as shown in figure 1. 




FiGUKE 1. Contour of integration of Hankel functions j H^l\z). 
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If \=u, llio part of the integral for H^^^ (z) from 
to 1 becomes, 

On the semicircle from +1 to — 1, if X=exp (id), 
this part is, 

1 f 

exp (iz sin d) exp {—vid)dd. (16) 

^ Jo 

From— 1 to— oo ^ if x^u'^e'"", this part is, 

^ r ^""p B (''"«)] '''"' '""p (-'"^^)'^^ (17) 

and together, 

Hl,^\z)=- Qxp (iz sin e—vid)de 
^ Jo 

+u'~'^ exp (— z/TT'Ol^/t^ (1<S) 
for Ke {z) > 0. 

AsT^^^) (z) = [Hi'^ (z^)]*' [Erdelyi, 1953], 

1 f 
^.^"^ (^)=- exp (— i;^ sin d+vid)dd 

TT Jo 

-1 — exp ^(^^~~) 0/"""^+^^""^ exp vTri)du 

(19) 

for Ke (2) > 0. Note that this method is applicable 
not only for v=\ and |, but also for complex v. 

The integrals in (18) and (19) may be evaluated by 
Gaussian quadrature ]Kopal, 1955]. In this method 
a finite integral is expressed as a sum: 



/ 

J a 



b M 

F(x) dx=J^ W.^F{yJ+e {M) (20) 



where e (M) is an error term which may, in general, 
be made arbitrarily small by increasing M where 
m=\, 2, 3, . . . M. 



ym=^^[(h—a)Xm+{h^ra)]. 



(21) 



The XnS are the Gaussian abscissas and M determines 
the number of values of x to be used in the quad- 
rature. The Gaussian weights and abscissas may 
be determined from the following: 



/- 



1 M 

f(x)dx=^ Hmj{Xm) 

-1 



TO=1 



1 



W„^=-{h—a)H„,. 



(22) 
(23) 



3 The * denotes the conjugate of the quantity. 



The xjs are the roots of the Legendre polynomials 
defined by 

{x'-l)'^=2'^7n\P,r,{x) 



dx'^' 

Po(x) = l 
Pi(x) = x 



(24) 



P^(x)=^x' 



3 



P,ix)-- 



35 
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Polynomials of higher degree are determined by use 
of the recursion formula : 

(m+l)P,,+i(^) + mP,,,_i(x) = (2m+l)xP^(x). (25) 

Upon determination of the roots, the weight co- 
efficients, H,ri, of the corresponding quadrature 
formulas arc evaluated as follows : 



Hm 



(1 XffJ [PM{Xjn)\ 



(26) 



Forty-eight weights and abscissas were used in this 
quadrature [Davis and Rabinowitz, 1956] approxi- 
mating tlie integrand in this case with a polynomial 
of ninety-fifth degree [Kopal, 1955]. 

Equations (18 and 19) are valid for Re(s)>0, 
but the values of the Hankel functions in the second 
and tliird quadrants may be found from the 
equations : 

Hi^\z)=[mi\z*)r, (27) 



sm VTT 



— e 



-...^J}UI!J^H?^{z), (28) 



and 



sm vw 



(29) 



^(2) (,,...) sin (1+m).. ^ 
sm vir 

, ^. sin mvT -r-m^ / X 

J^^viri _^ HI^\Z) 

sm VTT 

with m= — 1 in this instance [Watson, 1948]. 
3. Method of Iteration 



An iterative method with the above procedure for 
finding Hankel functions may be used in (^11) to 
find the roots ot (5). The iterative method described 
by Muller [1956], although it was developed to solve 
polynomial equations, may be used for any function 
which is analytic in the neighborhood of the roots. 
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The functions for r are from (11) : 

F{r) 
and 









3 


1 


|5.|>1. 
(31) 



Alternate expressions for F{t) are used so that d^ 
does not completely dominate the function and 
prevent convergence. 

Each iteration in the process of finding one root 
is obtained from the calculation of the nearer root 
which passes through the last three points of the 
function, /(r). The quadratic equation is, in 
general, complex; i.e., it possesses complex roots and 
complex coefficients. The function /(r) has the 
same ordinates as F{t) for n-\-l distinct abscissas 
as Tn, Tn-i . . . tq. Thc Lagrange interpolation 
formula througli the three points 7^-2, F{ri-2), 
Ti_i, F{Ti-i), and n, Fin) of F{t) is, 



/(r)=^2T^+^ir+A 



(32) 



where ^2, -4i, and Aq are evaluated from the re- 
quirement : 

/(r,_2) = F(r,_2), /(r,_0=i^(r.-i) and J{r)=F(u). 
Upon solving for ^2, ^1, ^0 and using the quantities 

Pi Pi-1 

and 8i=l+\i the Lagrange interpolation formula 
becomes a quadratic in X: 

J(r)=\^-^[F(u^2)\'i-F(u^,)\A 

+FiT,)X,]+Xd-'[F(T,^2)\'i-F(r,.^)8l 

+F(t,)(\,+8,)]+F(u), (33) 

Ti^i is found for the condition, /(r) = 0, solving for X 
and employing the relationship, 



X — Xi+i — " 



Ti-~Ti-l 



(34) 



Rationalizing the numerator of the standard quad- 
ratic formula, and solving for X<+i, 

^._ -2Fiu)8t 

9i±^g'i-^F{u)5MFiri-2)h-Firi-i)St + F{r,)] 

(35) 
where, 

g,=F(T,.,)\',-(FT,.,)8l+F(r,)(h+8,). (36) 

The sign choice in the denominator of X^+i (35) is 
resolved by selecting the value for the larger denomi- 
nator. This choice of X^+i gives r^+i the root of the 
LaGrange interpolation formula (32) which is nearer 
Tif Tof Ti, and 72, the initial value of r are estimated 



in this instance from the limiting roots Ts,o(^e=^) and 
Ts,ay{^e=°°) for any desired s [Johler et al., 1956], 
with the corresponding F(t) evaluated from (30 or 
31). 

The actual values of Ts,o and Ts,a, are used for s 
<4 while approximations for other s are determined 
from the following [Miller, 1946]: 






A.1 



-21/3 



48y! S&yU 
7 , 35 1 



where, 



mi 



'288yU 



2/,=^ (4s+3) 
y,=^ (4s+l). 



(37) 



(38) 



(39) 



(40) 



The iterative process is terminated with r^, the 
desired root when, 






<( 



(41) 



where e is a predetermined number. 

With procedures developed for finding Hankel 
functions and r's, Er may be calculated from (1) 
and (3). 

4. Discussion of Results 

In figures 2 and 4, describing effects of land and 
sea water, respectively, the transmitter on the sur- 
face, the receiver at varying heights above the sur- 

c 
face in wavelengths, X=:7 and for any distances along 

the surface of the earth, the ratio of the field aloft 

relative to that at the surface {h=0), j-,/ ' » ap- 

JlL {co, a) 

preaches a value of one as the height of the receiver 
approaches zero. For increasing heights of the re- 
ceiver, this ratio first exhibits a minimum less than 
one, although slight in some instances, and then in- 
creases in an exponential manner. The analytical 
behavior of the height-gain function has been dis- 
cussed in considerable detail by Wait [1956]. 

Characteristics of the amplitude ratio described 
above, however, with some differences, are present 
if both transmitter and receiver are elevated (fig. 6). 
For a transmitter height of 1 wavelength and vary- 
ing receiver height, for any distance along the sur- 
face of the earth, the amplitude ratio approaches a 
value less than 1 as the receiver height approaches 
zero. For a transmitter of height 10 wavelengths 
and varying receiver height, the amplitude ratio ap- 
proaches values greater than one which varies with 
the distance along the surface as the receiver height 
approaches zero. The amplitude ratios for these 
elevated transmitters increase very rapidly as ex- 
pected because there are two height gain factors 
increasing in an exponential manner. 

For the phase relationships (figs. 3, 5, and 7) cor- 
responding to the amplitude ratios (figs. 2, 4, and 6) 
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FiGUEE 2. Amplitude of the vertical electric field, Er, relative 
to surface value at various distances over land, with either 
transmitter or receiver elevated. 

In this figure as weU as those which follow hi and hz can be interchanged 



10^ 



10 



^ "^ 10" 



10" 



10-^ 



^ 


I I I ITM] II I' I IIHT] I I M Mill ' ' ' 

VERTICAL ELECTRIC DIPOLE SOURCE 

t =1110:!)^ sec 


rnri] 
1000-^ 


T IIII lllj 

200 mi 
/5OO / 


^-1000 mi - 
500 
200 . 


LEAD 






LAG 




/ 


/loookc^ / 


: 






100 kc/si 


a = 75 


/ 


_ 


.2=15 
a -- 0.005 
f^^= 47r(IO-'') 


1000 mi 100^// / 
500---^ X 1 / 
200-—^^ ' / 

'°°'^/^ 10000 kcA 7 

)Okc/%v'^ / -1000 mi 
/ / 500 
/ / 200 
/ ^/ 100 

/ -1000 kc/% 


"^1000 mi 
500 
200 
100 


_ 


IC 

/ 

/ 
/ 


r 


/ 








-^ 


" 




1 






' 



ALTITUDE,^ , WAVELENGTHS 
A 



Figure 3. Phase of the vertical electric field, Er, relative to 
surface value at various distances over land with either trans- 
mitter or receiver elevated. 
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Figure 4. Amplitude of the vertical elect? ic field, E^ relative 
to surface value at various distances over sea water, with 
either transmitter or receiver elevated. 
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Figure 5. Phase of the vertical electric field, Er, relative to 
surface value at various distances over sea water with either 
transmitter or receiver elevated. 



respectively, the difference of the phase aloft and 
that at the surface {h=0), t {cj,d,h) — t{o),d)<C<CO 
with a very slight minimum at low increasing heights. 
To illustrate details of the phase in this region, a 
logarithmic scale was used with t {o),d,h) — t {o),d)<^0 
graphed as lead with dashed lines, and t (o),d,h) 
—t{o:,d)y>0 designated as lag with solid lines (figs. 
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Figure 6. Amplitude of the vertical electric field, Er, relative 
to surface value at various distances over land, with both 
transmitter and receiver elevated. 

3, 5). The phase lag may be ambiguous by 27r 
radians in any case, which in microseconds is 10, 1, 
or 0.1 for 106, 1,000. and 10,000 kc/s, respectively. 
In figure 7, the lower solid lines represent lead for 
the transmitter at 1 wavelength, while the upper 
solid lines represent lag for the same case. For 10- 
wavelength height of transmitter, the dotted lines 
represent lead. 

5. Conclusions 

All of the amplitude ratio curves have a minimum., 
even though slight in some cases, and then rise ex- 
ponentially. For transmitter and receiver both 
elevated this ratio is very pronounced. The phase 
curves also exhibit a minimum., sometimes very 
slight, with t{o),d,h) — t(cA:,d) changing from a nega- 
tive to a positive quantity or from a phase lead to 
a phase lag relative to the surface. 

In applying the methods used in this paper, it may 
be concluded that the iterative method described is 
applicable for roots of Riccati's differential eq(5) to 
any desired accuracy consistent with the capacity of 
the electronic computer used. Also, the method for 
evaluating Hankel functions, as it is applicable for 
real or complex order, may be used for the many 
other applications of these functions. 

In particular, the height gain functions and the 
ground wave field can be evaluated at most any dis- 
tance or altitude of transmitter and/or receiver with 
the aid of the techniques presented, the ultimate 
limitation being the computer capacity and speed, 
and the failure of the approximation in (3). 
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Figure 7. Phase of the vertical electric field, Er, relative to 
surface value at various distances over land with both trans- 
mitter and receiver elevated. 
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